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Closing  Developments  in  Aerodynamic  Simulation 
with  Disjoint  Patched  Meshes 

Charles  K.  Lombard 
PEDA  Corporation 

1.  Introduction 

Recent  and  projected  developments  in  supercomputers,  numerical  grid  gen¬ 
eration  techniques  and  computational  algorithms  for  the  compressible  Euler  and 
Navier-Stokes  equations  portend  a  major  revolution  in  the  manner,  pace,  and  cost 
of  design  and  the  resulting  performance  of  aerodynamic  systems.  To  realize  these 
potential  benefits,  certain  closing  developments  in  computational  technique  must 
be  made  in  order  to  effect  a  highly  accurate,  reliable,  efficient  and  productive  sim¬ 
ulation  environment  for  aerodynamic  design  analysis. 

A  primary  need  of  the  developments  is  to  achieve  the  capability  for  a  user  to 
easily  and  rapidly  perform  flowfield  calculations  accurately  among  problems  of  dis¬ 
parate  and  realistically  complex  geometries.  The  natural  approach  to  realizing  this 
objective  with  comparatively  straightforward  extensions  of  existing  computational 
technology  is  through  the  use  of  systems  of  quadrilateral  patched  meshes. 

Such  systems  can  be  either/both  composite  (joined)  or  overset  (disjoint).  In 
the  former  case  adjacent  patches  share  a  common  boundary,  or  at  least  parallel 
boundaries  in  the  case  of  mesh  patch  overlap  for  purposes  of  applying  numerical 
boundary  conditions.  With  composite  meshes,  patch  boundaries  are  piecewise  fitted 
to  segments  of  physical  or  computational  boundaries  or  embedded  flow  structures. 
As  shown  by  Lombard,  et  al1,  composite  mesh  systems,  that  may  have  numerically 
useful  properties  of  geometric  continuity  across  patch  boundaries,  admit  topolog¬ 
ically  singular  global  meshes  that  have  the  capability  to  connect  computational 
regions  of  great  (really  any)  geometric  complexity.  However,  situations  exist  where 
multiple  mesh  topologies,  each  naturally  related  to  some  different  piece  of  geome- 
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try  or  flow  structure,  offer  greater  flexibility  and  accuracy  than  composite  meshing 
alone.  Examples  involve  multiple  bodies  that  may  have  relative  motion  or  weak 
shocks  that  bear  little  geometric  relationship  to  boundaries  of  the  flow  domain.  In 
such  cases  systems  of  disparately  oriented  and  at  least  partially  overset  grids  as 
proposed  by  Berger  and  Oliger2  allow  arbitrarily  high  resolution  of  all  features  of 
the  flowfield. 

To  make  efficient,  productive  use  of  patched  meshing  strategies  requires  a  body 
of  new  computational  tools  and  methodology  that  are  the  objectives  of  the  present 
research.  The  needed  factors  are:  (l)  a  simple  procedure  for  generating  patched 
computational  meshes  with  freedom  of  point  and  gradient  specification  on  all  patch 
boundaries;  (2)  improved  upwind  algorithm/numerical  boundary  condition  proce¬ 
dures  for  semi-autonomous  implicit  but  unconditionally  stable  conservative  coupling 
of  solutions  on  a  system  of  multiple  patched  meshes;  and  (3)  computer  graphics, 
particularly  a  simple  algorithm  for  constructing  contour  plots  on  systems  of  over¬ 
set  patched  meshes.  The  glue  that  is  to  tie  these  tools  together  in  a  simulation 
requiring  minimum  human  intervention  to  adapt  to  new  configurations  is  a  flexible 
parameter  controlled  multiple  mesh  data  structure.  An  important  objective  of  the 
program  is  to  test  the  evolving  techniques  in  appropriate  problems. 

2.  Research  Accomplishments 

In  the  first  year  of  an  anticipated  three  year  program  the  emphasis  was  placed 
on  the  most  crucial  and  challenging  objectives  -  patched  grid  generation  and  robust 
upwind  algorithm/boundary  procedures  for  rapid  relaxation  on  multiple  meshes. 

Algebraic  Grid  Generation 

The  concept  of  patched  meshing  in  which  complex  domains  are  broken  up  into 
many  geometrically  regular  and  topologically  rectangular  subdomains  leads  natu¬ 
rally  to  the  use  of  efficient  algebraic  techniques  for  the  construction  of  the  individual 
mesh  patches.  To  obtain  the  desired  smoothness  properties  over  the  global  mesh 
in  the  vicinity  of  patch  boundaries,  a  technique  that  permits  specification  of  point 
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distribution  and  gradient  on  all  boundaries  wets  devised.  The  technique  -  termed 
generalized  transfinite  interpolation3  -  makes  use  of  a  parameterized  general  cubic 
polynomial  for  the  coordinate  curves.  Regularity  of  the  mesh  is  obtained  by  em¬ 
ploying  continuous  distributions  of  the  parameters  of  the  curves  within  judiciously 
chosen  bounds  based  on  analysis.  Stretching  functions  such  as  that  of  Vinokur4  are 
used  to  distribute  points  and  blending  functions  are  used  to  distribute  parameters 
of  the  curves  between  lateral  boundaries. 

A  novel  feature  of  the  technique  is  the  introduction  of  the  corner  singularity 
from  analysis  to  govern  distribution  of  points  and  parameters  in  the  vicinity  of 
boundary  slope  singularities.  At  such  points,  the  method  thus  obtains  the  desired 
properties  of  mesh  smoothness  to  the  interior.  A  global  mesh  solution  obtained  with 
the  method  for  a  backward  step  problem  is  shown  in  Figure  1.  Here  the  solution 
was  generated  in  two  patches  one  containing  the  exterior  corner  and  the  other 
the  interior  corner.  The  solution  was  matched  analytically  at  the  patch  interface. 
Additional  background,  details  of  the  technique  and  other  numerical  results  are 
contained  in  reference  3. 

Upwind  Implicit  Relaxation  Algorithm/Boundary  Procedures 

Under  the  contract  we  have  devised  a  new  single  level  operationally  explicit 
but  effectively  implicit  algorithm  for  gasdynamics.  The  algorithm  is  particularly 
appropriate  for  multiple  patch  mesh  systems  because  each  solution  sweep  operation 
on  any  patch  is  decoupled  from  any  other.  Thus  the  method  is  not  only  very  storage 
efficient  and  simple  to  program  including  coupling  at  patch  boundaries,  but  also  can 
make  excellent  use  of  parallel  computing  in  several  straightforward  ways. 

Previously  the  Beam-Warming  factored  implicit  algorithm5  with  the  Baldwin- 
Lomax  thin  layer  viscous  approximation6  has  provided  the  basis  for  two  similar 
space  marching  (PNS)  procedures7,8  for  the  compressible  Navier-Stokes  equations. 
These  PNS  methods  which  are  highly  efficient  -  requiring  half  the  data  storage 
and  a  small  fraction  of  the  computer  time  of  two  level  time  dependent  methods  - 
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have  proven  effective  for  flows9,10  with  favorable  streamwise  pressure  gradient  or 
with  relatively  small  adverse  pressure  gradients.  However,  in  the  presence  of  strong 
adverse  pressure  gradient  such  as  occurs  in  a  wing  or  fin  root  regions  the  contem¬ 
porary  PNS  methods  suffer  numerical  stability  problems  and  may  infer  streamwise 
separation  even  where  separation  doesn’t  occur11.  In  such  unseparated  (perhaps 
weakly  separated)  regions,  numerical  stability  has  been  maintained  at  the  price  of 
employing  large  amounts  of  artificial  viscosity  with  a  resulting  loss  in  predictive 
accuracy  and  knowledge  of  the  actual  state  of  the  flow.  Where  strong  streamwise 
separation  occurs  the  methods  are  unstable  and  cannot  proceed.  Particularly  for 
the  increasingly  relevant  laminar  flow  situation  that  will  be  encountered  at  very 
high  altitude  by  aerodynamic  systems  such  as  orbital  transfer  vehicles  (OTV’s)  and 
space  shuttle,  streamwise  separation  becomes  a  likely  occurrence12  in  compression 
corners  associated  with  canopies,  pods,  flared  bodies,  wing  or  fin  roots  and  de¬ 
flected  control  surfaces.  Thus  a  more  general  technique  is  needed  that  is  inherently 
stable  for  all  types  of  upstream  influence.  At  a  minimum  the  mixed  elliptic  hyper¬ 
bolic  problem  requires  global  iteration,  preferably  with  type  dependent  differencing. 
Various  steps  in  this  direction  have  recently  been  taken  by  Rakich13,  by  Rubin  and 
Reddy14  and  by  Rizk  and  Chaussee15. 

Rakich13  utilized  global  iteration  with  Vigneron’s  Mach  number  dependent 
upwind  and  downwind  splitting8  of  the  pressure  gradient  for  the  streamwise  mo¬ 
mentum  equation.  The  approach  provides  an  improvement  in  both  accuracy  and 
stability  for  both  weakly  interactive  boundary  layer  flow  and  strongly  interactive 
flow  without  (significant)  streamwise  separation.  Global  iteration  is  carried  out  by 
marching  the  PNS  equation  repeatedly  only  along  the  downstream  direction. 

For  incompressible  flow  Rubin  and  Reddy14  (also  Lin  and  Rubin16  for  super¬ 
sonic  flow)  introduced  type  dependent  (upwind)  differencing  of  the  streamwise  ve¬ 
locity  along  with  pressure  splitting  in  an  implicit  method  involving  a  staggard  grid 
-  dependent  variable  location  scheme.  This  method  was  extended  to  compressible 
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flow  by  Khosla  and  Lai17.  The  method  admits  strong  interaction  with  streamwise 
separation  but  is  not  homogeneous  across  the  spectrum  of  Mach  number  regimes  of 
interest  here.  Other  PNS  related  methods  for  subsonic  flow  and  boundary  layer  flow 
with  streamwise  reversal  have  recently  been  reviewed  by  Brown18  in  conjunction 
with  a  new  staggered  grid  scheme. 

In  an  attempted  generalization,  Rizk  and  Chaussee15  presented  a  hybrid  tech¬ 
nique  of  space  marching  in  supersonic  zones  and  relaxation  with  a  Beam- Warming 
time  dependent  central  difference  method  over  zones  of  strong  upstream  influence 
with  streamwise  separation.  They  presented  two  variations  of  the  time  dependent 
method:  one  fully  implicit  requiring  two  levels  of  storage  and  matrix  inversion  pro¬ 
cedures  in  all  space  directions,  the  other  explicit  in  all  space  directions  but  the  thin 
layer  direction  and  requiring  only  one  level  of  storage  as  in  a  space  marching  algo¬ 
rithm.  Both  procedures  require  the  same  several  hundred  iterations  minimum  to 
reasonably  converge,  with  the  semi  explicit  procedure  requiring  substantially  less 
machine  time.  However,  both  procedures  are  inherently  unstable,  except  for  the 
use  of  artificial  dissipation.  The  marginal  stability  coupled  with  the  lack  of  type 
dependent  differencing  and  well  posed  boundary  approximations  all  contribute  to 
slow  convergence.  In  balance,  while  workable,  the  hybrid  approach  leaves  something 
to  be  desired  from  the  points  of  view  of  convenience,  computer  time  and  internal 
consistency  of  the  global  solution  procedure. 

A  new  globally  iterated  scheme  related  to  that  which  will  be  presented  here,  has 
been  presented  by  Moretti19.  In  his  procedure  for  steady  inviscid  flows,  the  Euler 
equations  are  cast  in  Riemann  variables  and  the  resulting  uncoupled  equations  are 
solved  by  integrating  each  Riemann  variable  separately,  sweeping  back  and  forth 
alternatively  along  each  coordinate  line.  The  coupling  between  the  equations  and, 
thus,  the  non-linearity  are  introduced  only  through  the  boundaries  and  the  updating 
of  state  after  complete  sweeps.  As  can  be  inferred  from  the  results  to  be  presented 
here,  the  reduced  coupling  inherently  limits  the  rate  of  convergence  of  Moretti’s 


method  and  the  procedure  is  not  extendable  to  the  compressible  Navier-Stokes 
equations. 

New  Universal  Single  Level  Scheme  CSCM-S 

The  CSCM  flux  difference  eigenvector  split  upwind  implicit  method20,21,22  for 
the  inviscid  terms  of  the  compressible  Navier-Stokes  equations  provides  the  natural 
basis  for  an  unconditionally  stable  space  marching  technique  through  regions  of 
subsonic  and  streamwise  separated  flow.  In  such  regions  the  split  method  can  be 
likened  to  stable  marching  of  each  scalar  characteristic  wave  system  in  the  direction 
of  its  associated  eigenvalue  (simple  wave  velocity).  In  supersonic  flow,  where  all 
eigenvalues  have  the  same  sign,  the  method  automatically  becomes  similar  to  the 
referenced  PNS  techniques  based  on  the  Beam- Warming  factored  implicit  method 
with  the  Baldwin-Lomax  thin  layer  viscous  approximation. 

Compared  to  contemporary  central  difference  methods,  the  CSCM  character¬ 
istics  based  upwind  difference  approximation  with  its  inherent  numerical  stability 
leads  to  greatly  reduced  oscillation  and  greater  accuracy  in  the  presence  of  captured 
discontinuities  such  as  shocks,  contacts  and  physical  or  computational  boundaries. 
The  correct  mathematical  domains  of  dependence  that  correspond  with  physical 
directions  of  wave  propagation  are  coupled  with  well  posed  characteristic  boundary 
approximations21  naturally  consistent  with  the  interior  point  scheme.  The  result 
is  faster  sorting  out  of  transient  disturbances  and  substantially  more  rapid  conver¬ 
gence  to  the  steady  state.  The  splitting  and  the  associated  time  dependent  implicit 
method  have  been  described  in  detail  in  references  (20)  and  (22)  for  quasi  1-D  and 
2-D  planar  or  axisymmetric  flow. 

In  the  following,  we  will  sketch  the  differences  between  the  time  dependent 
method  and  the  new  space  marching  technique  which  we  designate  CSCM-S.  The 
discussion  will  begin  with  the  quasi  1-D  inviscid  formulation,  present  some  results 
elucidating  the  properties  and  performance  of  the  method,  then  give  additional 
details  entering  into  multidimensional  inviscid  and  thin  layer  viscous  procedures 


and,  lastly,  present  early  2-D  solutions  obtained  with  the  new  single  level  scheme 
in  problems  solved  previously22  with  the  time  dependent  method. 

Quasi  1-D  Formulation 


The  general  jth  interior  point  difference  equations  for  the  time  dependent 
CSCM  upwind  implicit  method  is  written 

(/+  A+V  +  A-A)%  =  -A+Aq)  ”  -  A~ Aq)n  (l) 

3  -  1  3 

where  V  and  A  are  backward  and  forward  spatial  difference  operators.  In  the 
notation  the  interval  averaged  matrices  between  node  points  j  and  j  +  1  are  labeled 
j.  The  right  hand  side  of  equation  (l)  is  written  for  the  first  order  method.  Higher 
order  methods  in  space  are  given  with  results  in  references  20  and  22. 

Central  to  its  accurate  shock  capturing  capability,  the  CSCM  conservative  flux 
difference  splitting  has  the  “property  U”  put  forth  by  Roe23 
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(A+  +  A~)Aq)j  =  A F)j  =  F:  +  1  -  Fs  (2) 

Here  q  is  the  conservative  dependent  variable  vector  and  F  is  the  associated  flux 

/V 

vector.  The  matrices  A+  and  A~  are  the  splittings  of  the  CSCM  interval  averaged 
Jacobian  matrix  according  to  the  signs  of  the  averaged  eigenvalues.  Thus  in  the 
equation  for  the  jth  grid  point,  A+A^g)y_i  represents  stable  characteristic  spatial 
differencing  backward  for  positive  eigenvalue  contributions  and  A~A(q)j,  forward 
for  negative  ones. 

With  Sq  =  qn+1  —  qn ,  equation  (1)  defines  a  two  level  linearized  coupled  block 
matrix  implicit  scheme  that  can  be  solved  by  a  block  tridiagonal  procedure.  In 
reference  (22)  a  new  (DDADI)  approximately  factored  alternating  sweep  bidiagonal 
solution  procedure  for  equation  (1)  is  presented  that  is  shown  to  be  very  robust  and 
is  effectively  explicit,  i.e.  requires  only  a  decoupled  sequence  of  local  block  matrix 
inversions  rather  than  the  solution  of  the  coupled  set.  For  the  forward  sweep  the 
bidiagonal  solution  procedure  can  be  written 

(/  +  A+  -  A~)6q*j  =  RHS  +  A+6q*j_1 
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For  the  linear  problem,  i.e.  constant  coefficient  case  of  stability  analysis,  equation 
(3)  is  equivalent  to  the  single  level  space  marching  procedure 

{I  +  A+  -  A~)6q*,  =  A+q*-_1  -  A+qn  -  A~ Aq)H  (4) 

3  3 

Nonlinearity  enters  in  the  single  level  space  marching  form  (4)  in  that  at  each  step  of 
the  forward  sweep  the  matrices  A+  are  averaged  between  q* y_x  and  g”:  rather  than 
homogeneously  at  the  old  iteration  level  n.  Similarly,  a  companion  backward  space 
marching  sweep  that  is  symmetric  to  equation  (4)  and  that  is  intimately  related  to 
the  backward  sweep  of  the  alternating  bidiagonal  algorithm  of  reference  (22)  is 

(I  +  A+-  A~)6qj  =  -1+AO/-!  +  A~q)  -  A~ q^  1  (5) 

The  method  given  by  equations  (4)  and  (5)  is  von  Neumann  unconditionally  sta¬ 
ble  for  the  scalar  wave  equation.  The  analysis  shows  the  significance  of  DDADI 
approximate  factorization  in  rendering  both  the  forward  and  backward  sweeps  sep¬ 
arately  stable  regardless  of  eigenvalue  sign.  Consequently  as  the  local  Courant 
number  becomes  very  large,  the  robust  method  becomes  a  very  effective  (symmetric 
Gauss-Seidel)  relaxation  scheme  for  the  steady  equations,  a  fact  which  substantially 
contributes  to  the  very  fast  performance  that  will  be  demonstrated. 

At  a  right  computational  boundary  on  the  forward  sweep  we  solve  the  charac¬ 
teristic  boundary  point  approximation22 

(A+  +  A+)S,^=A\]J_i-A+,£  (6) 

9^+1  =  q  ^  and  at  a  left,  on  the  backward  sweep 

[A~  -  A~)6qi  =  A~qn^- A~qn  +  l  (7) 

Following  the  solution  of  equations  (6)  and  (7)  the  conservative  state  vector  is 
iteratively  corrected21  to  maintain  the  accuracy  of  prescribed  boundary  conditions 
while  not  disrupting  the  representation  of  the  computed  characteristic  variables 
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running  to  the  boundary  from  the  interior.  Analysis  of  a  model  system  with  upwind 
differenced  scalar  equations  and  coupled  boundary  conditions  was  related  to  the 
linearized  bidiagonal  scheme22  by  Oliger  and  Lombard25;  the  analysis  also  strongly 
supports  the  numerically  confirmed  robust  stability  of  the  present  nonlinear  method 
for  gasdynamics.  In  contract  to  reference  22,  a  useful  result  of  reference  25  is  that 
on  the  forward  sweep  there  is  no  need  for  a  predictor  step  at  the  left  boundary 
J  =  1;  to  the  solution  sweep  begins  at  J  =  2.  Similarly,  the  backward  sweep  begins 
at  J  =  N  —  1. 

With  the  updating  at  each  step,  where  in  equation  (5)  6qj  —  q™+1  —  <7y,  it  is 
clear  that  the  symmetric  pair  of  equations  (4)  and  (5)  serve  to  advance  the  solution 
two  pseudo  time  (iteration)  levels;  whereas,  the  linear  alternating  bidiagonal  sweep 
algorithm  of  reference  (22)  advances  the  solution  only  one  level.  To  maintain  con¬ 
servation  to  a  very  high  degree,  in  single  sweep  marching  in  supersonic  zones  we 
iterate  (at  least)  once  locally  at  each  space  marching  step.  The  local  iteration  serves 
to  make  the  eigenvectors  in  the  coefficient  matrices  consistent  with  the  advanced 
state  and,  thus,  provides  improved  accuracy  for  the  nonlinear  system.  It  appears 
effective  to  do  this  inner  iteration  everywhere,  i.e.  in  both  subsonic  and  supersonic 
regions,  as  the  number  of  global  iteration  steps  to  convergence  with  two  inner  iter¬ 
ations  has  been  found  reduced  by  a  factor  of  three  to  four.  Since  the  computational 
work  per  two  steps  is  about  the  same  for  the  single  level  and  two  level  schemes  and 
beyond  the  fact  that  one  saves  a  level  of  storage  in  the  space  marching  algorithm, 
the  question  arises:  Can  one  get  solutions  in  less  computational  work  through  faster 
convergence  with  the  nonlinear  space  marching  algorithm? 

One  Dimensional  Results 

First,  we  present  results  for  supersonic  flow  with  no  shock  in  Shubin’s  diverging 
nozzle.  In  purely  supersonic  zones,  the  experience  with  the  present  method  is  that 
the  solution  can  be  marched  accurately  in  one  global  iteration,  as  ought  to  be  the 
case.  Figure  2  shows  the  exact  solution  (in  solid  line)  and  the  computed  result 
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from  the  first  forward  sweep.  It  is  evident  that  the  method  correctly  predicts 
the  solution  to  plotting  accuracy  in  one  forward  sweep.  With  subsequent  sweeps 
the  error  (the  difference  between  the  exact  and  the  computed  solution)  reduces  to 
machine  accuracy  in  less  than  three  global  iterations.  In  fact,  by  increasing  (from 
two)  the  number  of  inner  iterations  on  the  solution  procedure  at  each  space  marching 
step,  convergence  to  prescribed  accuracy  can  be  guaranteed  in  one  forward  sweep. 
This  is  also  true  of  cor  temporary  locally  linearized  unsplit  methods  in  supersonic 
flow. 

With  the  globally  iterative  nonlinear  space  marching  formulation,  early  expe¬ 
rience  in  two  quasi  1-D  nozzle  problems  with  mixed  supersonic-subsonic  zones  is 
that  solutions  are  obtained  in  roughly  an  order  of  magnitude  fewer  iteration  steps 
than  had  been  required  with  the  previously  fast  pseudo  time  dependent  technique 
and  block  tridiagonal  solving. 

The  two  nozzle  problems  which  are  described  and  solved  by  Yee,  Beam  and 
Warming26  and  solved  with  the  CSCM  time  dependent  technique  in  references  (20) 
and  (21)  are  Shubin’s  diverging  nozzle  flow  and  Blottner’s  converging-diverging 
nozzle  flow.  Both  problems  involve  unmatched  overpressures  at  the  outflow  which 
result  in  internal  shock  terminated  supersonic  zones  and  subsonic  outflow.  For 
the  experiments  involving  flow  of  mixed  type  the  same  initial  data  given  by  Yee, 
Beam  and  Warming  -  a  linear  interpolation  between  inflow  and  outflow  values  for 
effectively  exact  solutions  of  the  problems  -  is  used  that  was  used  previously  with 
the  time  dependent  approach. 

For  flows  of  mixed  type,  in  Figures  3  and  4  respectively,  results  are  shown  for 
successive  forward  and  backward  sweeps  for  five  global  iteration  steps  with  Shubin’s 
and  Blottner’s  nozzle  flows.  In  both  cases,  the  exact  solution  as  given  by  Yee,  Beam 
and  Warming  is  shown  in  solid  line  and  the  present  computational  results  solved 
on  a  51  point  mesh,  in  boxes.  Blottner’s  nozzle  flow  is  shown  converged  after  10 
global  iteration  steps.  There  is  substantial  evidence  in  other  results  not  shown  that 


with  further  work  the  number  of  global  iterations  required  to  compute  flows  such 
as  Blottner’s  can  be  reduced  by  a  factor  of  two,  to  about  five. 

In  Figure  5,  we  show  a  subcritical,  i.e.  completely  subsonic,  flow  solution 
computed  in  only  two  global  iteration  steps  for  the  Blottner  nozzle  geometry  with 
different  inflow  conditions.  Here  the  exact  analytical  solution  derived  by  Venkata- 
pathy  is  shown  in  solid  line  and  our  computed  results  in  boxes. 

The  alternating  direction  sweeps  in  our  method  have  been  derived  directjy  out 
of  theory  for  solving  the  implicit  set  of  difference  equations.  However,  one  can  see 
mechanistically,  numerically  speaking,  that  omitting  the  backward  sweep  from  the 
pair  and  globally  iterating  only  with  the  forward  sweep  equation  (4)  will  result  in 
permitting  the  influence  of  a  subsonic  outflow  boundary  (or  interior  disturbance)  to 
propagate  upstream  only  one  grid  point  per  global  iteration.  In  such  a  case,  which 
relates  to  other  global  iteration  methods  found  in  the  literature  and  that  also  sweep 
only  in  the  main  flow  direction,  the  rate  of  convergence  is  greatly  inhibited  relative 
to  symmetric  sweeping  by  a  factor  of  order  roughly  the  number  of  grid  points  in  the 
subsonic  zone.  Mathematically,  this  inhibition  is  the  result  of  the  failure  to  include 
in  the  implicit  process  the  effect  of  the  eigenvectors  governing  upstream  influence 
but  to  treat  these  waves  explicitly  with  effective  CFL  unity. 

In  Figure  6  we  illustrate  the  progress  of  the  transient  solution  to  the  subcritical 
nozzle  problem  after  15  forward  sweeps,  with  the  backward  sweeps  omitted.  One 
can  clearly  see  that  the  wave  influence  of  the  outflow  boundary  has  progressed  only 
15  mesh  points  forward  of  the  outflow  boundary.  In  Figure  7  the  transient  solution 
is  shown  after  60  steps  which  is  beyond  one  characteristic  transit  time  (equivalent 
to  50  mesh  intervals)  for  the  upwind  wave  to  reach  the  inflow  boundary.  In  Figure 
8  we  show  the  history  of  the  RMS  error  in  the  primitive  variables.  The  solution  is 
found  to  converge  to  roughly  the  same  RMS  error  after  three  characteristic  times 
(150  steps)  as  the  solution  obtained  with  the  symmetric  alternating  sweep  sequence 
after  only  3  global  iterations. 
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Blottner’s  supercritical  nozzle  problem  which  involves  subsonic  inflow  acceler¬ 
ating  through  a  sonic  point  to  a  supersonic  zone  terminated  by  a  shock  to  subsonic 
outflow  is  the  most  computationally  demanding  of  the  test  cases  and  indicates  the 
capability  for  the  method  to  compute  simply  and  consistently  over  the  subsonic 
forebody  and  base  regions  of  blunt  bodies  in  supersonic  flow.  Thus  the  need  for 
separate  time  dependent  codes  will  be  obviated  by  this  new  method. 

Finally,  in  Figures  9  and  10,  we  present  the  convergence  history  for  the  present 
nonlinear  scheme  and  the  linearized  time  dependent  scheme  for  completely  subsonic 
and  supersonic  nozzle  flows.  The  x-axis  shows  the  number  of  iterations  each  scheme 
requires  to  reduce  the  exact  error  to  five  orders  of  magnitude  for  various  Courant 
numbers.  It  is  evident  that  the  present  scheme  converges  extremely  fast  at  all  CFL 
numbers  compared  with  the  method  based  on  the  linearized  block  tridiagonal  solver. 

Two  Dimensional  Formulation 

For  two  dimensional  flow,  assuming  a  marching  coordinate  £,  inviscid  terms 

B+V^  +  B-  A„  (8a) 

and 

Anq)k-i  —  B~Ar)q)k  (86) 

are  added  to  the  left  and  right  hand  sides  respectively  of  both  the  forward  and  back¬ 
ward  sweep  equations  (4)  and  (5).  For  viscous  flow,  second  centrally  differenced, 
thin  layer  viscous  terms  are  also  added  in  the  rj  direction  as  is  conventionally  prac¬ 
ticed,  e.g.  Steger27.  With  the  terms  for  the  r\  cross  marching  coordinate  direction, 
the  technique  now  becomes  an  implicit  method  of  lines.  Along  each  r\  coordinate 
line,  one  can  solve  the  equations  coupled  with  a  block  tridiagonal  procedure.  Alter¬ 
natively,  a  further  DDADI  bidiagonal  approximate  factorization  can  be  employed 
in  the  rj  direction  and  solved  either  linearly  as  in  reference  (22)  or  nonlinearly  as 
here  in  the  £  direction.  As  shown  in  the  quasi  2-D  numerical  experiments  of  ref¬ 
erence  (22),  DDADI  bidiagonal  approximate  factorization  is  stable  for  viscous  as 


well  as  inviscid  terms.  Finally  in  reference  (22)  there  is  a  relevant  discussion  of  the 
reduced  approximate  factorization  error  that  attends  using  DDADI  in  one  or  more 
space  directions.  A  variety  of  multidimensional  solution  strategies  derived  from 
DDADI  bidiagonal  approximate  factorization  in  a  simple  symbolic  algebra  based 
on  the  difference  stencil  are  presented  in  reference  24. 

While  extensions  to  3-D  are  not  given  in  detail  here,  we  note  such  extensions  to 
a  developing  two  level  3-D  CSCM  upwind  scheme  are  possible  and  will  be  adopted  in 
the  future.  For  3-D,  the  inviscid  terms  similar  to  (8a)  and  (8b)  will  again  be  added 
for  the  f  cross  flow  direction.  For  efficiency  in  solving  each  resulting  marching  plane, 
the  implicit  operators  for  the  r\  and  f  coordinate  directions  can  be  approximately 
factored  using  DDADI  as  described  in  reference  (22)  for  two  space  directions. 

Two  Dimensional  Results 

We  present  results  for  a  45°  -  15°  axisymmetric  transonic  nozzle  flow  previously 
studied  experimentally  by  Cuffel,  Back  and  Massier28  and  computationally  by  Cline, 
Prozan,  Serra  and  Shelton  (all  referenced  in  (28))  and  ourselves22.  In  Figure  11  we 
show  results  after  10  steps  of  an  early  computation  run  at  a  local  CFL  number  of 
20  with  the  present  first  order  single  level  scheme.  Except  for  the  addition  of  an 
error  correction  procedure21  to  counter  numerical  inflow  boundary  condition  drift, 
a  factor  which  has  improved  the  present  solution  in  the  vicinity  of  the  axis,  the 
effectively  converged  results  found  here  are  the  same  as  those  given  for  the  two 
level  scheme  in  reference  22.  (As  long  as  the  problem  has  a  unique  solution,  the  two 
schemes  must  give  equivalent  results  since  the  right  hand  side  difference  equation 
sets,  including  boundary  approximations,  are  the  same.) 

For  the  solution  given  in  Figure  11,  we  noted  a  very  rapid  rate  of  reduction 
in  residual,  three  orders  of  magnitude  in  ten  steps.  This  compares  with  60  steps 
given  in  reference  (22)  for  the  solution  obtained  with  the  two  level  scheme.  The 
rapid  convergence  found  in  this  transonic  problem  for  the  CSCM-S  method  with 
viscous  terms  provides  the  reasonable  expectation  of  similar  fast  results  to  be  ob- 


tained  without  viscous  effects.  Thus  the  method  in  multidimensions  appears  to  have 
attractive  potential  for  an  improved  transonic  Euler  solver  as  well  as  Navier-Stokes 
solver. 

Next,  we  present  first  order  inviscid  and  viscous  results  for  an  inlet  problem 
shown  in  Figure  12.  The  pressure  contours  for  the  first  order  inviscid  solutions  are 
shown  in  Figure  13.  Figure  14  shows  the  first  order  viscous  results.  The  viscous 
computation  shows  the  presence  of  the  leading  edge  shock.  The  flow  structure  com¬ 
pares  very  well  with  the  theoretical  (for  the  inviscid  case)  and  other  computational 
results.  In  Figure  15,  the  inviscid  and  viscous  wall  pressure  are  compared  with 
the  exact  solution  (inviscid).  Figure  16  shows  the  convergence  history  of  the  RMS 
residue  of  all  the  conservative  flow  variables  for  the  inviscid  problem  solved  at  CFL 
=  100  with  4  inner  iterations  at  every  axial  location.  For  the  inviscid  case,  only 
forward  marching  was  carried  out  and  backward  marching  was  omitted.  The  solu¬ 
tion  has  converged  for  practical  purposes  at  the  end  of  the  first  sweep.  The  residue 
reaches  machine  accuracy  in  10  iterations.  In  a  later  paper29,  we  show  the  residual 
reduction  versus  inner  iteration  number  in  single  sweep  solutions  for  supersonic  flow 
and  compare  results  with  contemporary  PNS  procedures. 

Patched  Mesh  Boundary  Procedures 

In  previous  work  Lombard,  et  al21,22  and  Oliger  and  Lombard25  gave  stable 
implicit  procedures  for  computing  the  solution  at  external  boundaries  of  a  com¬ 
putational  domain.  Those  procedures  generalized  the  work  of  Kenzer  to  matrix 
coupled  linearized  boundary  conditions  complementing  a  set  of  advective  difference 
equations  (associated  with  well  posed  characteristics)  to  the  boundary. 

Under  the  present  contract  we  have  explored  the  problem  of  implicitly  cou¬ 
pling  at  interior  patch  boundaries  the  global  solution  on  a  system  of  multiple  patch 
meshes.  The  approach  we  have  taken  in  this  research  is  numerical  experimenta¬ 
tion  among  a  number  of  boundary  treatment  approximations  to  the  solution  of  the 
continuous  domain  problem.  For  comparison  purposes  with  previous  single  mesh 


results,  numerical  experiments  were  performed  with  the  well  tested  two  level  pseudo 
time  relaxation  CSCM  scheme  at  the  interior  points  of  the  mesh  patches. 

Our  first  experiments,  to  be  described  here,  dealt  with  breaking  a  single  co¬ 
ordinate  line  into  segments  and  solving  sequentially  on  each  the  equations  for  a 
quasi  2-D  viscous  compressible  flow.  The  model  problem,  with  which  we  experi¬ 
mented  in  reference  22  for  an  uninterrupted  mesh,  is  a  transient  pipe  flow  resulting 
from  an  initially  nonequilibrated  pressure  between  the  axis  and  wall  boundaries. 
Three  kinds  of  cases  were  run  with  the  two-level  linearized  implicit  procedure;  all 
featured  sequential  solving  on  patches  with  at  least  one  point  of  overlap.  Case  1 
had  frozen  boundary  data,  obtained  from  the  solution  on  neighboring  meshes  in  an 
operationally  explicit  manner.  Specifically  at  left  and  right  first  computed  interior 
points  of  a  patch,  we  solved  the  bidiagonal  equations  respectively 

{I  +  A+-  A~A)6q2  =  A+q *  -  (A+  -  A~)q ”  -  A~q ”  (9a) 

and 

(I-A-  +  A+V)SqN^  =  A+qN1  2-(A+-  A~)qNn_  J  -  A-qnN  (96) 

Here  the  symbols  n,  n  +  1  indicate  the  “frozen”  boundary  data  to  come  from  the 
solution  on  the  interior  of  an  adjoining  and  partially  overlying  patch  may  be  either 
at  the  old  or  new  iteration  level  in  the  global  procedure.  In  Case  1  the  solution  on 
all  the  grids  was  effectively  updated  at  the  same  time. 

Case  2  featured  reverse  sequential  cycling  of  the  two  level  time  dependent  solu¬ 
tion  procedure  through  the  grids  on  alternate  global  steps.  In  the  forward  sequence 
the  right  interior  patch  boundaries  were  (a)  frozen  or  (b)  computed  with  one-sided 
characteristic  boundary  conditions  obviating  any  change  in  the  characteristic  data 
from  outside  (to  the  right  of)  a  patch.  The  left  interior  patch  boundaries  inherited 
implicit  data  ( 6q )  from  the  solution  on  the  computed  patch  to  the  left.  In  the 
backward  sequence  all  the  roles  were  reversed  including  the  directions  of  forward 
elimination  and  back  substitution  in  the  tridiagonal  matrix  inversion  procedure. 


Case  3  featured  cycling  through  the  patches  in  the  predictor  (forward  elimi¬ 
nation)  step  of  the  solution  procedure,  inheriting  implicit  left  boundary  data  as  in 
Case  2.  Then  the  patches  were  cycled  through  in  reverse  order  on  the  corrector 
(back  substitution)  step.  The  result  (save  for  interpolation  if  data  points  of  the 
grids  were  interlaced)  is  identically  equivalent  to  solving  on  an  uninterrupted  single 
patch  mesh. 

The  results  of  the  three  sets  of  experiments  were  comparable  for  local  time  steps 
based  on  constant  CFL  number  up  to  about  50.  Beyond  that  the  rate  of  convergence 
of  the  solutions  for  Cases  1  and  2  diminished  and  at  sufficiently  high  CFL  number 
failed  to  converge.  The  effort  involved  in  Case  2  with  the  two  level  scheme  was  not 
rewarded  relative  to  the  simple  procedure  of  Case  1.  Within  the  framework  of  the 
single-level  symmetric  Gauss-Seidel  implicit  relaxation  scheme,  however,  Case  2a  is 
operationally  no  more  difficult  than  Case  1  and  more  closely  approximates  Case  3. 
Case  2b  has  a  consistency  problem  that  inhibits  firm  convergence. 

Figure  17  compares  rms  residual  (density)  convergence  history  for  Cases  1  and 
3  with  three  mesh  segments  with  two  cells  of  overlap  and  run  at  CFL  25.  As  one 
might  expect,  the  effectively  uninterrupted  mesh  procedure  is  found  to  converge 
(two  to  three  times)  faster  for  about  the  first  50  steps  through  the  major  transient. 
After  that  the  performance  is  comparable.  The  performance  of  the  frozen  boundary 
treatment  Case  1  with  five  segments  at  CFL  25  was  found  to  be  no4  materially  worse 
than  with  three  segments,  but  the  performance  degradation  with  increasing  CFL 
was  found  to  proceed  faster  with  increasing  number  of  patches. 

From  the  results  of  the  quasi  2-D  experiments  and  extrapolation  from  expe¬ 
rience  with  the  single  level  scheme,  we  observe  that  the  simple  boundary  proce¬ 
dure  with  frozen  conservative  variable  data  taken  from  the  solutions  on  adjacent 
meshes  and  coupled  implicitly  by  alternating  direction  sequential  solving  through 
the  patches  has  robust  stability  to  sufficiently  high  CFL  number  to  yield  a  rate  of 
convergence  meeting  our  needs. 
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4.  Interactions 

The  research  described  in  this  report  has  to  this  point  been  partially  presented 
in  the  form  of  a  paper  on  algebraic  grid  generation  by  Vinokur  and  Lombard  (refer¬ 
ence  3),  a  SIAM  meeting  paper  (reference  25)  by  Oliger  and  Lombard  on  boundary 
procedures  for  bidiagonal  alternating  sweep  schemes  and  a  Computational  Fluid 
Dynamics  Seminar  at  NAS  A- Ames  Research  Center  by  Lombard  on  the  Univer¬ 
sal  Single  Level  Implicit  Algorithm.  The  latter  invoked  tremendous  interest  and 
discussion. 

The  research  is  about  to  spawn  three  other  papers  on  the  single  level  relaxation 
algorithm  -  references  24  and  29  and  the  unreferenced  paper 

Lombard,  C.K.,  Venkatapathy,  E.  and  Bardina,  J.:  Forebody  and  Baseflow  of 
a  Dragbrake  OTV  by  an  Extremely  Fast  Single  Level  Implicit  Algorithm,  AIAA- 
84-1699,  Snowmass,  Colorado,  June  1984. 

5.  New  Discoveries 

The  Universal  Single  Level  Algorithm  for  the  compressible  Euler  or  Navier- 
Stokes  equations  is  a  new  discovery  in  numerical  methods  that  promises  to  result 
in  substantial  efficiences  in  data  storage,  programming,  machine  time  and  human 
productivity. 
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Figure  1.  Portion  of  an  algebraically  generated 
computational  mesh  for  a  flow  domain  containing  an 
external  and  an  internal  corner. 


Figure  2.  Shubin's  diverging  nozzle  supersonic 
flow  solution  developed  in  one  forward  sweep  from 
supersonic  initial  data. 
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Figure  3.  Shubin's  diverging  no;zle  flow  solution  developed  in  alternating  forward 
and  backward  sweeps,  one  each  per  global  iteration  step  for  five  steps, 
square  -  computed  results  ,  line  -  exact  solution 
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Figure  3.  continued 


continued 
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Figure  4.  Blottner's  converging-diverging  supercritical  nozzle  flow  solution 
developed  in  alternating  direction  sweeps  for  ten  global  iteration 
steps.  square  -  computed  results  ,  line  -  exact  solution 


Figure  4.  continued 


Figure  4.  continued 


Figure  6.  Solution  to  subcritical  nozzle  problem  after 
15  global  iterations  with  forward  marching  sweeps  only. 
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Figure  If.  Convergence  history  of  the  RMS  of 
the  residuals  for  the  inviscid  first  order  inlet 
problem  with  4  inner  iterations  per  global  sweep. 
Note  at  the  end  of  the  first  sweep  the  residual 
has  reduced  and  the  solution  has  converged  for  all 
Dractical  purposes. 


Figure  17.  Test  on  three  patches  of  implicit 
stability  and  rate  of  convergence  of  (circles) 
computed  boundary  point  operator  differencing 
to  frozen  data  in  adjacent  patches;  (solid 
line)  effectively  continuous  grid  method. 


